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Abstract 

Large  amplitude,  time-domain,  wave-body  interactions 
are  studied  in  this  paper  for  problems  with  forward 
speed.  Both  two-dimensional  strip  theory  and  three- 
dimensional  computation  methods  are  shown  and  com- 
pared by  a number  of  numerical  simulations.  In  the 
present  approach,  an  exact  body  boundary  condition 
and  linearized  free  surface  boundary  conditions  are 
used.  By  distributing  desingularized  sources  above  the 
calm  water  surface  and  using  constant-strength  flat 
panels  on  the  exact  body  surface,  the  boundary  in- 
tegral equations  are  solved  numerically  at  each  time 
step.  Once  the  fluid  velocities  on  the  free  surface 
are  computed,  the  free  surface  elevation  and  poten- 
tial are  updated  by  integrating  the  free  surface  bound- 
ary conditions.  After  each  time  step,  the  body  surface 
and  free  surface  are  regrided  due  to  the  instantaneous 
changing  wetted  body  geometry.  The  desingularized 
source  method  applied  on  the  free  surface  produces 
non-singular  kernels  in  the  integral  equations  by  mov- 
ing the  fundamental  singularities  a small  distance  out- 
side of  the  fluid  domain.  Constant-strength  flat  panels 
are  used  to  model  bodies  with  any  arbitrary  shape. 

Extensive  results  are  presented  to  validate  the  ef- 
ficiency of  the  present  methods.  These  results  include 
the  added  mass  and  damping  computations  for  a mod- 
ified Wigley  hull  and  a S-175  hull  with  forward  speed 
using  both  two-dimensional  and  three-dimensional  ap- 
proaches. Heave  force  and  pitch  moment  time  his- 
tories of  a S-175  hull  due  to  a heave  motion  using 
a time-domain  body-exact  strip  theory  are  presented. 
Diffraction  forces  acting  on  a modified  Wigley  hull  due 
to  a linear  head  sea  incoming  wave  using  fully  three- 
dimensional  method  are  also  obtained.  All  the  com- 
putational results  are  compared  with  experiments  or 
other  numerical  solutions. 

1 INTRODUCTION 

The  accurate  prediction  of  wave-induced  motions  and 
loads  are  very  important  in  ship  and  offshore  design, 
which  requires  knowledge  of  the  maximum  value  of  var- 
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ions  design  parameters  such  as  ship  motion  amplitude, 
bending  moments,  impact  forces,  etc. 

A mixed  Euler-Lagrange  time-stepping  scheme 
(MEL)  was  first  introduced  by  Longuet-Higgins  & 
Cokelet  (1976)  for  solving  two-dimensional  fully  non- 
linear water  wave  problems.  Since  then,  MEL  methods 
have  also  been  successfully  used  to  solve  fully  nonlin- 
ear, three-dimensional  wave  and  wave-body  interaction 
problems  (Dommermuth  & Yue,1987;  Xue  et  ah,  2001; 
Cao,  Schultz  & Beck, 1990,  1991;  Cao,  1991;  Cao,  Lee 
fe  Beck,  1992;  Scorpio  et  ah,  1996).  The  problems 
with  MEL  computations  are  the  instabilities  of  the  free 
surface  and  wave  breaking.  The  instabilities  can  of- 
ten be  eliminated  by  improved  numerical  techniques, 
but  wave  breaking  is  a natural  phenomenon  that  is 
expected  to  occur  in  any  large  body  motion  or  wave 
situation.  Computations  normally  are  forced  to  stop 
when  wave  breaking  occurs.  Various  techniques  have 
been  proposed  to  continue  the  computations  after  wave 
breaking,  but  they  are  still  not  robust  and  can  lead  to 
nonphysical  solutions. 

A compromise  between  fully  nonlinear  computa- 
tions and  linear  theory  is  the  so-called  body-exact  ap- 
proach. In  the  body-exact  problem,  the  body  bound- 
ary condition  is  satisfied  on  the  instantaneous  wetted 
surface  of  the  body  while  the  linearized  free  surface 
boundary  conditions  are  retained.  In  order  to  solve  for 
the  hydrodynamic  forces  due  to  large  body  motions 
in  the  body-exact  problem,  a time-domain  approach 
is  preferred.  A method  to  deal  with  the  exact  body 
boundary  condition  using  a time-domain  free  surface 
Green  function  has  been  developed  by  Beck  & Magee 
(1990)  for  a submerged  body  performing  arbitrary  mo- 
tions. Other  researchers  such  as  Lin  & Yue  (1990) 
also  have  successfully  obtained  results  for  a surface- 
piercing body  problem  using  the  time-domain  free  sur- 
face Green  function  method. 

In  order  to  develop  computationally  fast  seakeep- 
ing calculations,  but  still  retaining  the  important  non- 
linearities,  a time-domain  body-exact  method  using 
Rankine  sources  has  been  developed.  Earlier  work 
on  two-dimensional  body-exact  computations  has  been 
given  by  Zhang  and  Beck  (2006,  2007).  They  solved 
two-dimensional  large  amplitude  radiation  and  diffrac- 
tion problems  including  water  entry  and  exit.  Com- 


parisons  with  the  other  numerical  calculations  and  ex- 
periments were  good. 

In  this  paper,  we  have  continued  to  develop  two- 
dimensional  and  three-dimensional,  body-exact  com- 
putational models  together  with  using  a Rankine 
source  Green  function  to  compute  the  large  amplitude 
body  motions  in  waves.  The  exact  body  boundary  con- 
dition is  applied  on  the  instantaneous  wetted  hull  sur- 
face under  ~ = 0.  For  the  large  amplitude  compu- 

tations, the  instantaneous  submerged  body  surface  is 
repanelized  at  each  time  step  to  capture  the  nonlinear- 
ities due  to  the  changing  geometry.  The  linearized  free 
surface  boundary  conditions  are  used. 

To  compute  the  three-dimensional  problem  re- 
sults, two  approaches  have  been  followed.  The  first 
is  a strip  theory  in  which  a series  of  two-dimensional 
problems  at  different  cross-sections  are  solved  at  each 
time  step.  Once  the  boundary  value  problems  have 
been  solved  at  each  section,  Radial  Basis  Functions 
(RBF)  are  used  to  compute  the  longitudinal  deriva- 
tives of  the  velocity  potential  on  the  body  surface.  In 
the  second  approach,  the  three-dimensional  problem 
is  solved  completely  and  the  derivative  of  the  veloc- 
ity potential  on  body  is  computed  directly.  For  the 
large  amplitude  radiation  problem  with  forward  speed, 
both  methods  are  used  to  compute  the  hydrodynamic 
forces,  and  both  sets  of  results  are  compared  with  ex- 
periments. 

As  will  be  formulated  in  the  following  sections,  the 
three-dimensional  solver  includes  the  forward  speed  ef- 
fects both  in  the  free  surface  boundary  conditions  and 
the  pressure  calculations.  The  derivatives  of  the  poten- 
tial on  the  hull  can  be  easily  computed  once  the  source 
strengths  have  been  obtained  by  solving  the  boundary 
value  problem.  In  the  strip  theory  calculations,  the 
boundary  value  problem  is  solved  at  different  cross- 
sections  and  the  forward  speed  correction  is  only  in- 
cluded in  the  pressure  computations  by  computing  the 
derivative  of  potential  in  the  x-direction  on  the  hull 
using  Radial  Basis  Functions.  The  computation  time 
in  solving  the  three-dimensional  problem  is  relatively 
more  expensive  than  solving  the  problem  using  a time- 
domain  strip  theory.  The  accuracy  of  those  results  will 
be  shown  and  compared  in  the  following  sections. 

2 MATHEMATICAL  FORMULATION 
2.1  Three-Dimensional 

A boundary  value  problem  for  a vessel  travelling  in 
deep  water  is  solved.  The  vessel  moves  with  forward 
speed  U(t),  and  may  be  undergoing  unsteady  oscilla- 
tions in  its  six  degrees  of  freedom.  The  fluid  is  assumed 
to  be  ideal  and  the  flow  irrotational.  Two  coordinates 
systems  will  be  employed:  the  x0  system  is  fixed  in 
space,  and  x system  is  fixed  to  the  mean  position  of 
the  ship  and  is  moving  with  forward  speed  U(t)  along 
the  track  of  the  ship.  The  boundary  value  problem 
is  solved  in  the  right  hand  moving  coordinate  system 


Figure  1:  Problem  sketch  and  reference  frames 


( x,y,z ),  as  shown  in  Figure  1.  The  x-axis  points  out 
the  bow,  and  the  z-axis  upward.  The  origin  is  at  the 
calm  water  line  at  mid-ship. 

A velocity  potential  is  introduced  to  describe  the 
fluid  motion  by  using  the  above  assumptions  such  that 
the  fluid  velocity  can  be  expressed  as  the  gradient 
of  potential  function,  V{x,i)  = VI?  = V(— U(t)x  + 
6(x.  y,z,t)),  where  <j>  is  the  perturbation  velocity  po- 
tential. 

The  form  of  the  linearized  free  surface  boundary 
conditions  are: 


dr/ 

~di 

Ut 


-gri  + U ■ - — 


(1) 

(2) 


In  order  to  use  an  Euler-Lagrange  free  surface  time 
stepping  scheme,  the  free  surface  boundary  conditions 
can  lie  rewritten  in  a more  convenient  form: 


§ = ^ + U-Vri  + v-Vr,  (3) 

§ = —gg  + U-Vc/)  — — +v  - V(j)  (4) 

Ot  p 


where  jj  = gj  + v ■ V is  the  time  derivative  follow- 
ing a fluid  particle  along  a prescribed  path.  The  ve- 
locity of  the  particle  in  moving  coordinate  system  is 
v = (u,v,0).  Here,  u and  v are  the  prescribed  veloc- 
ities for  the  horizontal  plane  motion  of  the  free  sur- 
face collocation  points,  u is  allowed  to  move  with  the 
translation  velocity  —U,  and  v is  prescribed  so  that 
collocation  points  move  on  the  given  paths  around  the 
body. 

This  formulation  has  the  desired  effect  of  restrict- 
ing collocation  points  from  passing  through  the  body 
boundary.  An  alternative  free  surface  boundary  condi- 
tion would  be  to  set  v = (0,0,0),  using  equations  (1) 


and  (2).  In  this  case  ^}  = jjf  and  the  collocation  points 
are  fixed  on  the  calm  water  surface,  tracking  the  free 
surface  elevations  at  points  in  the  x — y plane  like  wave 
probes.  The  advantage  of  this  scheme  is  that  no  regrid- 
ing of  the  free  surface  after  each  time  step  is  required 
if  the  body  is  wall-sided.  The  disadvantage  is  that  the 
free  surface  slope,  Vry,  must  be  computed  numerically. 
V?7  may  be  difficult  to  compute  accurately  in  a three- 
dimensional  problem,  especially  at  the  boundaries  of 
the  free  surface  domain  where  one  sided  differencing 
must  be  used.  The  free  surface/body  intersection  line 
is  where  this  problem  would  be  most  detrimental. 

All  the  velocity  potentials  should  satisfy  the 
Laplace  equation  under  the  assumption  of  ideal  po- 
tential flow.  By  applying  Green’s  theorem  and  using  a 
Rankine  source  Green  function,  the  velocity  potential 
can  be  written  as 

<!>=[[  G(x,$)a(£)ds  (5) 

J J SpUSs 

where  G = — L_-  a js  the  source  strength  on  the 

r(x,f) 

boundary;  Sp  is  the  calm  water  free  surface;  and  Spit) 
is  the  instantaneous  wetted  body  surface.  The  exact 
body  boundary  condition  is 

n-V(t>=U0(t)ni(t)  + VH-n(t)  on  SB(t)  (6) 


a velocity  potential  z,  t).  In  the  fluid  domain,  if) 
satisfies  Laplace’s  equation 

VV  = 0 (10) 

On  the  mean  free  surface,  the  linearized  free  sur- 
face boundary  conditions  are  imposed 

Ct-if’z  = 0 on  z = 0 (11) 

ipt  + gC  = 0 on  z = 0 (12) 

where  z = C(vP)  is  the  free  surface  elevation,  g is  the 
acceleration  due  to  gravity.  On  the  instantaneous  body 
boundary,  no  normal  flux  is  permitted 

= VN  + U0tan(rj5)N3  - U0tan(r)6)N2  on  SB 

" (13) 

where  N is  the  two-dimensional  unit  normal  vector, 
and  is  positive  out  of  the  fluid.  Vjv  is  the  instanta- 
neous velocity  in  the  normal  direction  including  rota- 
tional effects.  The  last  two  terms  are  two-dimensional 
corrections  when  there  is  a non-zero  pitch  angle,  r/ 5, 
or  non-zero  yaw  angle,  jy,.  In  the  far  field,  a radiation 
boundary  condition  is  imposed  such  that  there  are  no 
incoming  waves;  also,  the  water  is  assumed  deep  and 
the  potential  vanishes  as  z — >■  —00.  The  initial  condi- 
tions at  t = 0 are 


where  U0(t ) is  the  time-dependent  translating  veloc- 
ity of  the  body  in  the  x direction;  n is  the  inward  unit 
normal  on  the  body  surface(out  of  fluid);  n 1 is  the  com- 
ponent of  unit  normal  in  the  x direction;  and  VB  is  the 
motion  velocity  including  rotational  modes  of  a point 
on  the  ship  surface. 

The  perturbation  potential  <j>  must  satisfy  the  ra- 
diation boundary  condition  such  that  there  must  lie 
no  incoming  waves  and,  in  the  deep  water  problem,  <p 
vanishes  as  z — > —00.  The  initial  conditions  at  t = 0 
can  be  written  as: 

<f>  = tj>t  = 0 in  the  fluid  domain  (7) 

After  solving  the  boundary  value  problem,  the 
force  and  moment  can  be  determined  by  using: 

P = j j pnds  (8) 


V;  = ij>t  = 0 in  the  fluid  domain  (14) 

At  each  time  step  a mixed  boundary  value  problem 
must  be  solved;  the  potential  is  given  on  the  free  surface 
and  the  normal  derivative  of  the  potential  is  known 
on  the  body  surface.  In  terms  of  the  desingularized 
sources  (as  described  in  the  next  section)  above  the  free 
surface  and  sources  distributed  on  the  body  surface, 
the  potential  at  any  point  in  the  fluid  domain  can  be 
given 


JL  r 

V(x)  = £>(&)  ln|x-&|  + / <7(0G(x;|)di 

i=1  sB 


(15) 


where  SB  represent  the  instantaneous  wetted  body  sur- 
face. |x  — £i|  represents  the  distance  between  any  point 
in  the  fluid  domain  and  the  desingularized  source  point. 
G(x:£)  is  a Rankine  source  Green  function 


M = 


p(f  x n)ds 


(9) 


G(x;  £)  = lnr  (16) 

r = |x-£|  (17) 


2.2  Two-Dimensional  Strip  Theory 

The  three-dimensional  seakeeping  problem  can  be  ap- 
proximated by  solving  a series  of  two-dimensional  prob- 
lems at  cross-sections  (stations)  of  the  vessel.  This  is 
valid  for  long  and  slender  ships.  The  two-dimensional 
approach  used  here  is  based  on  the  work  of  Zhang  and 
Beck  (2007).  The  method  follows  the  same  conventions 
as  described  for  the  three-dimensional  case.  The  two- 
dimensional  flow  at  each  station  can  be  described  by 


where  r is  the  distance  between  a source  point  and  a 
collocation  point;  £ is  the  source  point  on  the  body 
boundary. 

Once  the  source  strengths  are  found,  tjj  can  be 
evaluated  by  (15),  and  the  velocity  on  the  body  Viy 
can  been  obtained.  The  total  pressure  is  given  by 
Bernoulli’s  equation 

P=~P{1H+Im2-U°lte+9Z)  (18) 


Under  the  slender  body  assumption,  the  forces 
acting  on  the  body  can  be  approximated  by  integrating 
(18)  over  the  instantaneous  submerged  body  surface, 
which  can  be  written  as 


F = J dx  J pN  dl 

(19) 

J dx  J p(rxN)  dl 

(20) 

L 


3 NUMERICAL  METHODS 
3.1  Three-Dimensional 


< Tj  is  the  vector  of  unknown  source  strengths.  bj  is  the 
vector  of  boundary  conditions,  6,  = <j>0i  oil  the  free 
surface,  bj  = U0(t)ni(t)  + Vjj  ■ n(t ) on  the  wetted  body 
surface. 

Once  the  source  strengths  are  known,  the  particle 
velocities  on  the  free  surface  can  be  computed.  Then 
the  free  surface  conditions  are  updated  by  using  a 4th- 
order  Runge-Kutta  scheme.  Meanwhile,  the  velocity 
potential  on  the  body  surface  can  also  be  obtained. 
By  using  a central  time  differencing  scheme,  the  pres- 
sure acting  on  the  body  can  be  calculated  so  that  the 
force  can  be  obtained  by  integrating  the  pressure  over 
the  instantaneous  wetted  surface.  The  pressure  on  the 
body  surface  is  given  by  the  Bernoulli's  equation 


We  distribute  the  desingularized  sources  above  the 
calm  water  surface.  The  desingularized  distance  is  cal- 
culated according  to  the  formula  Ds  = LdS 2 (Lee, 
1992),  where  Ds  is  the  desingularized  distance,  S is  the 
local  grid  area,  and  Lci  = 1 is  a desingularized  parame- 
ter. Constant-strength  flat  panels  are  used  on  the  body 
surface.  The  instantaneous  submerged  body  surface  is 
discretized  into  N quadrilateral  elements  over  which 
the  source  strength  is  assumed  constant.  The  nonpla- 
nar  quadrilaterals  are  mapped  to  planar  elements  by 
fitting  the  corner  points  in  a least-squares  sense  (Hess 
and  Smith,  1964;  Newman,  1986).  A boundary  inte- 
gral equation  can  be  solved  for  the  unknown  panel  and 
isolated  source  strengths. 


p=  -?[f -vdi +k£>2  + (#)’  +(£)W] 


Td± 

dx 


dx 


dy 


' dz ' 


(23) 

The  dynamic  forces  acting  on  the  ship  in  mode  j due  to 
a unit  motion  in  mode  k can  be  found  by  integrating 
the  consequent  pressure  over  the  instantaneous  sub- 
merged body  surface. 


F; 


jk 


’ !LdSl 


9<j>k  rTd(t>k 

' dt  dx 


+ + (24) 


In  this  paper,  we  compute  the  derivatives  of  the  veloc- 
ity potential  on  the  hull  directly. 
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where 


3.2  Two-Dimensional  Strip  Theory 

3.2.1  Body  and  Free  Surface  Modeling 

Applying  the  boundary  conditions  and  using  equation 
(15),  the  following  integral  equations  can  be  solved  to 
determine  the  unknown  source  strength. 


(21) 


^V(&)ln|xc -&|+  f ff(£)G(xc,|)dZ  = V»(xc)  xc  € I> 

r _ J 

(25) 


rF 


rB 


Xi 

= (xj, yt,Zi ),  a collocation  point 

!>(«<)— 

b 

+ 

•u7> 

1 

u 6 

fii 

= (nXi,  nVi,  nZi),  unit  normal  pointing  into  liody 

Tf 

“ ■ - .j 

rB 

Sj 

= a panel  integration  surface  on  the  body 

M 

= number  of  panels  on  the  body  surface 

N 

= total  number  of  unknown  source  strengths 

where 

= a source  point 

9G(xc,Z) 

dn 


dl  = x(xc)  xc  6 
(26) 


When  the  above  summations  (21)  are  applied  at  N 
collocation  points  (a?,:),  an  N x N linear  system  results, 


AijOj  — bj 


(22) 


xc  = a point  on  the  real  boundary 
X = the  given  normal  velocity  on  the  body 
if}  = the  given  potential  on  the  free  surface 
I>  = the  free  surface 
Tg  = surfaces  on  which  x is  known 


where  the  influence  matrix  Ajj  = Gjj  = l/(|i?f  — |) 

for  collocation  points  on  the  free  surface  and  Ajj  = The  integral  mixed  boundary  value  equations,  (25) 

Hi  ■ for  collocation  points  on  solid  boundaries;  and  (26),  can  be  discretized  to  form  a system  of  linear 


equations.  On  the  free  surface,  desingularized  sources 
are  distributed  outside  of  the  domain  such  that  the  the 
source  points  never  coincide  with  collocation  or  node 
points,  avoiding  singularities.  Isolated  sources  are  used 
rather  than  a distribution,  reducing  the  complexity  of 
the  influence  matrix.  The  desingularized  distance  is 
given  by  the  square  root  of  the  local  mesh  size. 

The  free  surface  is  broken  up  into  an  inner  and 
outer  region.  The  inner  region  spans  four  wavelengths 
on  either  side  of  the  body,  where  the  wavelength  is 
given  by  the  dispersion  relation,  A = given  a fre- 
quency of  oscillation,  u>.  To  resolve  the  radiated  waves, 
30  nodes  are  distributed  per  wavelength  in  the  inner 
region.  The  outer  region  acts  as  a numerical  beach  to 
prevent  wave  reflection.  20  nodes  are  distributed  over 
80  wavelengths  with  exponentially  increasing  spacing, 
as  determined  by  Lee  (1992). 

The  body  surface  is  modeled  using  panels,  which 
are  more  suitable  for  arbitrarily  shaped  bodies.  The 
resulting  discretized  mixed  boundary  value  equations 
are  given  by 


Nf 
3 = 1 


hi  I 


ln|x£  -£f| dl  = ip(xB) 


(27) 
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Ldl  = x(xf ) 


(28) 


In  matrix  form,  equations  (27)  and  (28)  are  solved 
using  LU  decomposition.  Once  the  source  strengths 
are  known,  the  fluid  velocity  on  the  free  surface  and 
pressure  on  the  exact  body  can  be  computed. 


3.2.2  Time  Evolution 


The  body  is  repanelized  at  each  time  step  using  a “rub- 
ber band”  technique.  The  number  of  panels  remains 
the  same  for  each  station,  and  they  are  stretched  to 
where  the  exact  body  intersects  the  calm  water  sur- 
face, z = 0.  This  requires  a correction  in  the  dip/dt 
term  in  equation  (18),  such  that 


Ip*  — Ip*  At 
At 


)-vt  ■ ViP* 
(29) 

where  if b*  is  the  velocity  potential  at  time  step  f,  and 
v*  is  the  moving  node  velocity  due  to  repanelization. 

The  free  surface  elevation  and  potential  are  up- 
dated using  the  kinematic  (11)  and  dynamic  (12)  free 
surface  boundary  conditions.  Time  stepping  is  done 
using  a 3rd-order  Adams-Bashforth  scheme  as  shown 
below. 
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The  updated  free  surface  and  potentials  are  used 
to  start  the  mixed  boundary  value  problem  at  the  next 
time  step.  To  ensure  consistent  free  surface  resolution 
in  time,  the  free  surface  nodes  are  relocated  to  a distri- 
bution consistent  with  the  original  distribution  (based 
on  the  location  of  the  body  intersecting  the  free  sur- 
face). The  values  of  £ and  ip  are  interpolated  using 
cubic  splines  to  the  new  distribution. 

3.2.3  Section  Exit  and  Entry 

One  of  the  main  reasons  for  using  the  body-exact  ap- 
proach is  the  ability  to  get  accurate  solutions  for  large- 
amplitude  motions,  which  are  non-linear.  This  includes 
allowing  sections  to  come  out  of  and  go  back  into  the 
water.  This  is  done  simply  by  stopping  calculations  as 
a section  comes  out  of  the  water,  then  re-initializing 
ip  and  ipt  to  zero  at  that  section.  As  the  section  en- 
ters the  water,  the  two-dimensional  problem  is  solved 
as  usual.  The  velocity  potentials  are  used  in  pressure 
calculations  only  after  the  solution  steadies,  after  three 
time  steps. 

3.2.4  Forward  Speed  Corrections 

To  solve  the  full  three-dimensional  problem  with  for- 
ward speed,  classical  strip  theory  needs  a correction 
for  the  the  interaction  between  sections  along  the  lon- 
gitudinal direction.  With  non-zero  forward  speed,  the 
pressure  equation  (18),  including  the  moving  node  cor- 
rection, becomes 

P = -P(jf  + ||VV7|2  - Uo ^ „ • ViP  + gz ) (32) 

There  is  a necessity  to  find  the  dip/dx  term  in  the 
pressure  calculation  when  there  is  a non-zero  forward 
speed.  The  method  implemented  here  is  the  use  of 
radial  basis  functions.  This  can  only  be  done  once  the 
two-dimensional  problem  is  solved  at  each  section.  The 
velocity  potentials  on  the  body  must  be  known  at  each 
section.  The  details  of  this  method  are  discussed  in  the 
following  section. 

3.2.5  Radial  Basis  Functions 

The  forward  speed  problem  can  be  solved  by  using  a 
radial  basis  function  (RBF)  to  approximate  the  dij)/dx 
term  in  the  pressure  calculation.  Using  this  method  to 
numerically  solve  partial  differential  equations  is  dis- 
cussed by  Kansa  (1999).  A more  general  discussion  of 
RBFs  is  presented  by  Buhmann  (2000). 


Once  the  two-dimensional  problem  is  solved  for 
each  section,  the  velocity  potentials  at  nodes  on  the 
entire  body  are  given  as  a function  of  the  Euclidean 
distance  between  the  nodes. 


9. i(x)  = .9  ( 1 1 x — Xj||) 

(33) 
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where  g is  one  of  many  possible  basis  functions.  The 
final  equation  (35)  for  the  coefficients  oq  is  to  ensure 
uniqueness.  These  expansion  coefficients  are  found  by 
setting  up  a system  of  iV  + 1 equations,  such  that 


,9i(xi)  •••  9iv(xi) 

: gj(xi)  : 

9i(xjv)  •••  9n(xn) 
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and  the  N + 1 system  can  be  rewritten  in  matrix  form 


Ha  = tjj  (36) 

with  ^>(xjv+i)  = 0.  The  interpolation  expansion  coef- 
ficients can  be  found  after  matrix  inversion 


q = H~1i>  (37) 

These  coefficients  can  be  used  to  interpolate  the 
data  in  three-dimensional  space  and  to  find  derivatives 
at  any  point.  The  x-velocity  of  the  fluid  is  found  by 
taking  the  partial  derivative  of  the  basis  function  with 
respect  to  the  x-direction. 


dx 


N 

Y,aj9j(xi) 


3=1 


(38) 


Several  basis  functions  were  tested  and  compared 
before  implementation  in  the  code.  Of  the  many  func- 
tions available,  the  three  shown  below  exhibit  the  qual- 
ity of  a “sphere  of  influence",  such  that  points  closer 
to  the  node  of  interest  have  a larger  influence  on  that 
node 


9j(x)  = e ci  Exponential  Spline  (39) 


(40) 


Cj  as  a fraction  of  the  ship  length,  L 


Figure  2:  RBF  Cj  comparison 


c.  as  a fraction  of  the  ship  length,  L 


Figure  3:  RBF  c.j  comparison 
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g:j  ( x)  = (1+  l|X~2Xjl1  )-V2  Inverse  MQ  (41) 

Using  a defined  analytic  velocity  potential  as  a 
test,  the  candidate  functions  were  compared  in  ac- 
curacy of  interpolation  and  approximation  to  partial 
derivatives.  The  analytic  potential  was  of  the  form 
ip  oc  x2,  with  some  variation  for  nodes  on  the  same  sta- 
tion (x-location).  Another  function,  ip  <x  x2  + x3,  was 
tested  as  well.  The  approximate  partial  x-derivatives 
were  compared  to  the  known  solutions.  The  Gaussian 
equation  (40)  and  Inverse  Multi-Quadratic  equation 
(41)  proved  to  be  the  more  accurate  functions.  Their 
accuracy  depends  on  the  “sphere  of  influence”  term,  c,j. 
Several  tests  cases  were  run  to  determine  an  the  opti- 
mum value  of  this  variable.  They  involved  computing 
the  partial  x-derivative  of  the  RBF  approximation  and 
comparing  it  to  the  known  solution.  Both  the  Gaussian 
and  Inverse  Multi-Quadratic  functions  were  evaluated. 
Mean  error  and  standard  deviation  were  considered  as 
benchmarks  of  comparison.  The  percentages  (based  on 
the  mean  absolute  value  of  the  exact  solutions)  of  mean 
error  and  standard  deviation  are  shown  as  functions  of 
Cj  in  Figures  2 and  3.  They  indicate  that  using  approx- 
imately L/ 4 yields  good  results  for  both  functions. 

Implementing  this  radial  basis  in  the  code  resulted 
in  occasional  instabilities,  especially  for  extreme  am- 


Gaussian 


Figure  5:  Diffraction  force  (heave  force)  acting  on  a 
Wigley  hull  III,  Fn  = 0.3,  X/L  = 1.0,  /?  = 7 r 


Figure  4:  Desingularized  sources  distribution  above 
calm  water  and  flat  panels  discretization  on  the  hull 

plitude  motions  at  higher  frequency.  The  source  of  the 
instability  is  the  dip/dx  term.  It  is  a result  of  very 
high  condition  numbers  in  the  matrix  inversion  rou- 
tines, which  use  LU  decomposition.  Several  attempts 
were  made  to  correct  this,  but  some  instabilities  re- 
mained. The  Exponential  Spline  function,  equation 
(39),  although  less  accurate,  had  much  smaller  condi- 
tion numbers.  The  final  selection  for  use  in  the  code 
was  an  exponential  function  of  power  1.5,  as  shown  in 
equation  (42).  This  function  was  much  more  accurate 
than  the  Exponential  Spline  while  being  better  condi- 
tioned for  matrix  inversion. 

I lx~ x / 1 1 1-5 

,gi(x)  = e (42) 

The  results  were  again  verified  against  analytic  ve- 
locity potentials.  This  RBF  gives  accurate  approxima- 
tions to  di/j/dx,  while  remaining  stable  when  imple- 
mented in  the  real  code  under  real  conditions. 

4 RESULTS  AND  DISCUSSION 

Numerical  convergence  tests  for  both  two-dimensional 
and  three-dimensional  problems  have  been  shown  in 
previous  papers  (Zhang  and  Beck,  2007:  Zhang,  2007). 

4.1  Three-Dimensional  Results 

Figure  4 shows  the  desingularized  source  distribution 
above  the  z = 0 and  the  flat  panel  distributions  on  a 
modified  Wigley  hull.  Figure  5 and  Figure  6 show  the 
force  time  history  results  of  the  wave  diffraction  prob- 
lem for  a modified  Wigley  hull  at  Fn  = 0.3  in  head 
seas.  The  incident  wavelength  X/L  = 1.0,  where  L is 
the  ship  length  of  the  modified  Wigley  hull;  A is  the  in- 
cident wave  length:  (3  is  the  angle  of  wave  propagation 
measured  from  the  positive  sense  of  the  x-axis. 


Figure  6:  Diffraction  force  (pitch  moment)  acting  on  a 
Wigley  hull  III,  Fn  = 0.3,  X/L  = 1.0,  f3  = 7r 


4.2  Two-Dimensional  Strip  Theory  Results 

Figures  7 through  10  show  the  force  time  series  for 
forced  heave  motions  of  the  S-175  containership.  Posi- 
tive heave  displacement  indicates  the  vessel  is  above 
it’s  calm  water  position.  The  particulars  of  the  S- 
175  are  given  in  Table  1.  The  simulations  are  for 
small  and  large  amplitudes  of  heave,  at  low  and  high 
frequencies  of  oscillation.  The  simulations  were  run 
at  zero  and  varying  forward  speeds  corresponding  to 
Fn  = 0.0,0.15,0.3.  The  radiation  force  can  be  seen  at 
these  various  speeds. 

The  pitching  moment  is  most  sensitive  to  forward 
speed.  As  expected,  the  radiated  force  is  largest  for 
the  high  frequency  case.  Non-linear  behavior  can  be 
seen  in  the  large  amplitude  motion,  where  the  transom 
stern  section  of  the  body  exits  and  enters  the  water.  It 
is  especially  evident  in  the  high  frequency  case,  as  can 
be  seen  in  Figure  10. 


Lpp(m) 

175 

B(m) 

25.4 

D(m) 

15.4 

T(m) 

9.5 

A(t) 

24742 

Table  1:  Parameters  of  S-175  hull 


Heave  Force  due  to  Heave:  Amp=0.1  meters,  co=0.2654  rad/s 


0 20  40  60 

time  [sj 

Pitch  Moment  due  to  Heave:  Amp=0.1  meters,  ©=0.2654  rad/s 


0 20  40  60 

time  [s] 


Figure  7:  S-175  0.1  m forced  heave  at  0.2654  rad/s 


Heave  Force  due  to  Heave:  Amp=0.1  meters, co=1, 187  rad/s 


0 5 10  15 

time  [s] 

Pitch  Moment  due  to  Heave:  Amp=0.1  meters,  0=1.187  rad/s 


0 5 10  15 

time  [s] 


Figure  8:  S-175  0.1  m forced  heave  at  1.1870  rad/s 


Heave  Force  due  to  Heave:  Amp=4  meters,  ©=0.2654  rad/s 


0 20  40  60 

time  [s] 

Pitch  Moment  due  to  Heave:  Amp=4  meters,  ©=0.2654  rad/s 


0 20  40  60 

time  [s] 


Figure  9:  S-175  4.0  m forced  heave  at  0.2654  rad/s 


Heave  Force  due  to  Heave:  Amp=4  meters,  <0=1. 187  rad/s 


0 5 10  15 

time  [s] 

Pitch  Moment  due  to  Heave:  Amp=4  meters,  ©=1. 187  rad/s 


0 5 10  15 

time  [s] 


Figure  10:  S-175  4.0  m forced  heave  at  1.1870  rad/s 
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11:  Wigley  III  Coefficients,  Heave  due  to  Heave 
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Figure  13:  Wigley  HI  Coefficients,  Heave  due  to  Pitch 
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Figure  12:  Wigley  III  Coefficients,  Pitch  due  to  Heave  Figure  14:  Wigley  III  Coefficients,  Pitch  due  to  Pitch 


4.3  Comparisons 

Figures  11  through  14  display  the  heave  and  pitch 
added  mass  and  damping  coefficients  for  the  mod- 
ified Wigley  hull  at  Fn  = 0.3,  due  to  heave  and 
pitch  motions.  This  modified  Wigley  hull  corresponds 
to  Journee’s  ‘Wigley  III”,  as  referenced  in  Journee 
(1992).  Both  the  fully  three-dimensional  and  two- 
dimensional  strip  theory  codes  were  used  and  the  re- 
sults are  compared  the  experimental  results  of  Journee. 
Most  of  the  coefficients  compare  well.  The  first  excep- 
tion is  the  damping  in  pitch  due  to  heave,  £53  in  Figure 
12.  The  three-dimensional  results  are  high.  The  other 
exception  is  pitch  added  mass  due  to  pitch,  A55  in  Fig- 
ure 14.  The  current  methods’  values  are  low  compared 
to  experiments.  This  was  found  to  be  the  case  for  most 
computational  methods  when  comparing  to  Journee’s 
findings  for  these  coefficients. 

Figures  15  through  18  display  the  added  mass 
and  damping  coefficients  for  the  S-175  containership 
in  heave  and  pitch,  due  to  heave  and  pitch  motions. 
The  simulations  were  performed  with  forward  speed 
corresponding  to  Fn  = 0.275.  Both  the  fully  three- 
dimensional  and  two-dimensional  strip  theory  codes 
were  used  and  the  results  are  compared  with  SWAN 
version  2.2  as  obtained  by  personal  correspondence 
with  Prof.Sclavounos.  In  general,  the  coefficients  agree 
well.  Again,  an  exception  is  the  damping  in  pitch  due 
to  heave,  £53  in  Figure  16.  These  discrepancies  are 
consistent  with  the  modified  Wigley  hull  findings.  The 
two-dimensional  results  for  forced  pitch  motions  do  not 
agree  as  well  as  they  did  for  the  case  of  the  modified 
Wigley  hull.  Forward  speed  affects  pitch  motions  more- 
so  than  heave  motions,  and  the  S-175  shows  more  non- 
linearities  at  the  bow  and  stern. 

4.4  Conclusions 

Two-dimensional  and  three-dimensional  body-exact 
computational  models  are  formulated  and  developed 
in  this  paper.  Numerical  computations  of  the  hydro- 
dynamic  coefficients  for  a modified  Wigley  hull  and  a 
S7-175  hull  are  compared  with  experiments  and  other 
numerical  solutions.  The  comparison  shows  the  capa- 
bilities and  efficiency  of  the  present  two-dimensional 
and  three-dimensional  body-exact  models.  It  is  indi- 
cated by  comparisons  that  three-dimensional  free  sur- 
face effects  are  more  important  in  computing  the  hy- 
drodynamic coefficients  due  to  the  pitch  motion,  while 
time  domain  strip  theory  can  achieve  relatively  faster 
simulation. 

The  methods  discussed  here  are  similar  in  that 
they  model  the  ship’s  orientation  in  the  water  exactly 
and  utilize  the  time-domain  approach.  However,  they 
are  drastically  different  when  solving  interesting  prob- 
lems of  surge  and  forward  speed  (i.e.  calm  water  re- 
sistance). The  three-dimensional  approach  inherently 
captures  these  details,  while  the  two-dimensional  re- 
quires some  sort  of  forward  speed  correction  term.  The 


tradeoff  is  that  the  two-dimensional  strip  theory  is 
computationally  less  expensive,  typically  by  a factor 
of  8 ~ 10. 
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